REPORT  DOCUMENTATION 


Public  reoofting  ouraen  tor  ihis  coiieaion  ot  imormaiion  is  esnmaiefl  to  average  i  nour 
gathennq  ana  maintaining  the  data  needed,  ana  comoieting  ana  reviewing  me  collectiot 
colleaton  ot  intonmation.  including  suggestions  tor  reducing  tnis  ouraen,  to  Washington 
Davis  Highway,  Suite  1204.  Arlington.  VA  22202*4302.  and  to  me  Otiice  of  Manageme 


1.  AGENCY  USE  ONLY  (Leave  Wan/c) 


4.  TITLE  AND  SUBTITLE 


2.  REPORT  DATE 

January  1998 


AFRL-SR-BL-TR-98- 


wed 

')704-0188 

ig  existing  oata  sources. 

'  any  other  asoect  ot  this 
Reports.  1215  Jefferson 
)ton.  DC  20503. 


3.  REPORT  TYPE  AND  DATES  COVERED  ^ 

Technical,  9/96  -  12/97 


5.  FUNDING  NUMBERS 


Fracture  Phenomena  and  Interface  properties  of 
Composites:  A  Fundamental  Approach 


6.  AUTHOR(S) 

Arif  Masud  and  John  Botsis 


F49620-96-1-0467 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRE5S(ES) 

Dept,  of  Civil  and  Materials  Engineering 
University  of  Illinois  at  Chicago 
842  West  Taylor  Street 
Chicago,  IL  60607 


9.  SPONSORINGAVIONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AFOSR 

110  Duncan  Avenue,  Suite  B115 
Bolling  AFB  DC  20332-0001 


11.  SUPPLEMENTARY  NOTES 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSORINGAIAONITORING 
AGENCY  REPORT  NUMBER 


19980803  102 


12a.  DISTRIBUnON/AVAILABIUTY  STATEMENT 

12b.  DISTRIBUTION  CODE 

Approved  for  public  release } 
distribution  unlimited. 

13.  ABStBfKCT  (Maximum  200  words) 

This  report  presents  a  numerical  approach  to  study  the  effects  of  bridging  and  shielding  fibers 
on  reduction  in  Stress  Intensity  Factor  (SIF)  at  the  crack  tip.  The  variational  foundation 
of  the  formulation  is  based  on  a  modified  Hu-Washizu  variational  principle.  The  crack  is 
assumed  to  grow  in  a  self-similar  manner  and  an  energy  approach  is  undertaken  to  evaluate 
SIF  at  the  crack  tip.  A  perfect  bond  is  assumed  between  the  fiber  and  the  matrix  to  account 
for  the  interaction  of  the  two  phases.  It  is  observed  that  the  area  around  the  crack  tip  where 
the  fibers  strongly  influence  SIF  extends  to  about  three  fibers  ahead  of  the  crack  and  five 
fibers  in  the  wake  of  the  crack.  Numerical  simulations  also  show  that  there  is  a  certain  limit 
beyond  which  increasing  the  ratio  of  the  modulus  of  elasticity  of  the  fibers  w.r.t.  the  matrix 
may  not  help  in  reducing  the  SIF  at  the  crack  tip.  A  detailed  study  on  the  inter-phase 
stress  fields  around  the  reinforcing  fibers  shows  that  stronger  fibers  have  a  larger  zone  of 
influence  and  can  potentially  effect  the  stresses  around  the  neighboring  fibers.  A  study  on 
temperature  change  indicates  that  thermal  gradients  result  in  a  redistribution  of  the  stress 
field  and  thus  can  have  pronounced  effect  on  the  stress  field  at  the  crack  tip. 


14.  SUBJECT  TERMS 


Discrete  fiber-matrix  system,  finite  elements,  shielding 
effect,  bridging  effect,  thermal  effect 


17.  SECURITY  CLASSIRCAT10N  18.  SECURITY  CLASSIFICATION 
OF  REPORT  OF  THIS  PAGE 

Unclassified  Unclassified 


19.  SECURITY  CLASSinCATION 
OF  ABSTRACT 

Unclassified 


15.  NUMBER  OF  PAGES 


18.  PRICE  CODE 


20.  UMITATION  OF  ABSTRACT 

UL 


standard  Form  298  (Rev.  2-69) 
PfMCfibed  by  ANSI  Std.  Z39-18 


NSN  7540-01 -280-5500 


ABSTRACT 


Fibrous  and  laminated  composites  have  become  an  important  class  of  ma¬ 
terials  that  is  increasingly  being  used  in  engineering  structures.  However,  unlike 
the  common  engineering  materials,  characterizing  the  fracture  behavior  of  fibrous 
composites  is  a  difficult  task.  In  this  work  we  have  adopted  a  numerical  approach 
to  study  the  effects  of  bridging  and  shielding  fibers  on  reduction  in  Stress  Intensity 
Factor  (SIF)  at  the  crack  tip.  This  report  presents  a  finite  element  formulation  of 
elasticity  to  model  elastic  fracture  in  composites  with  long  aligned  fibers.  We  have 
employed  B-bar  type  approach  which  is  applicable  to  compressible  as  well  as  nearly 
incompressible  material  systems.  The  variational  foundation  of  the  formulation  is 
based  on  a  modified  Hu-Washizu  variational  principle.  The  crack  is  assumed  to 
grow  in  a  self-similar  manner  and  an  energy  approach  is  undertaken  to  evaluate 
the  stress  intensity  factor  at  the  crack  tip.  A  perfect  bond  is  assumed  between  the 
fiber  and  the  matrix  to  account  for  the  interaction  of  the  two  phases.  The  effects  of 
strong  intact  bridging  fibers,  fiber  ahead  of  the  crack  tip,  and  that  of  temperature 
variation  on  the  reduction  in  stress  intensity  factor  at  the  crack  tip  have  been  inves¬ 
tigated.  Various  numerical  results  are  presented  to  show  fiber-fiber  and  fiber-crack 
interaction. 
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SIGNIFICANT  FINDINGS 


We  have  adopted  a  novel  computational  approach  to  investigate  the  bridging 
and  the  shielding  mechanism  and  their  effect  on  the  stress  field  at  crack  tip  without 
any  other  dissipative  mechanisms.  The  following  main  points  can  be  drawn  from 
the  results  of  the  present  numerical  studies. 

1.  The  effect  of  the  bridging  fibers  on  SIF  is  much  larger  than  that  of  the 
fibers  ahead  of  the  crack  tip. 

2.  The  area  around  the  crack  tip  where  the  fibers  strongly  influence  the  stress 
intensity  factor  extends  to  about  three  fibers  ahead  of  the  crack  and  five  fibers  in 
the  wake  of  the  crack. 

3.  As  for  the  effect  of  the  modulus  of  elasticity  of  the  fibers  on  SIF,  the 
results  indicate  that  for  values  greater  than  Ej/Ejn  =  15  the  effect  is  relatively 
small.  Making  fibers  too  strong  with  respect  to  the  matrix  (beyond  a  certain  limit) 
may  not  help  in  reducing  the  stress  intensity  factor  at  the  crack  tip. 

4.  We  have  also  studied  the  inter-phase  stress  fields  around  the  reinforcing 
fibers  that  show  steep  gradients  normal  to  the  interface.  It  is  seen  that  stronger 
fibers  have  a  larger  zone  of  influence  and  can  potentially  effect  the  stresses  around 
the  neighboring  fibers.  The  present  study  also  confirms  that  homogenization  of  the 
composite  system  cannot  predict  the  discrete  stress  distribution  at  crack  tip  and 
fiber-matrix  interfaces,  and  thus  can  lead  to  erroneous  conclusions. 

5.  The  effects  of  temperature  change  have  also  been  studied  in  the  presence 
of  a  bridging  fiber.  These  preliminary  results  indicate  that  thermal  gradients  can 
result  in  a  redistribution  of  the  stress  field  and  thus  can  have  pronounced  effect  on 
the  stress  field  at  the  crack  tip.  It  also  shows  the  importance  of  a  coupled  thermo¬ 
mechanical  theory  for  possible  elasto-plastic  yielding  of  the  matrix  material  at  the 
tip  of  the  crack.  Such  elasto-plastic  deformation  would  considerably  alter  the  stress 
field  in  a  localized  zone  in  the  viscinity  of  the  crack. 
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1.  Introduction 


Composite  materials,  e.g.,  metal,  ceramic  and  polymer  reinforced  with  long 
aligned  fibers  have  shown  significant  promise  in  the  design  of  advanced  aerospace 
systems.  Their  high  strength  to  weight  ratio  and  the  fiexibility  to  tailor  make 
these  materials  with  performance  that  is  superior  to  their  constituent  materials 
are  the  main  attributes.  However,  unlike  the  common  engineering  materials,  the 
fracture  behavior  of  composites  has  not  been  well  understood.  Thus,  understanding 
strength  and  fracture  behavior  of  these  materials  has  been  the  subject  of  intense 
investigation  in  the  past  several  years.  Because  of  the  fact  that  the  microstructures 
of  real,  commercially  available  composite  materials  are  very  complex  and  often 
difficult  to  characterize,  it  has  been  difficult  to  reveal  basic  mechanisms  as  well 
as  elementary  events  that  lead  to  fracture  and  failure  in  these  advanced  material 
systems.  As  a  result,  the  lack  of  a  rigorous  failure  criteria  in  composites  made  of 
advanced  fibrous  materials  has  ebbed  away  the  confidence  in  using  these  materials 
in  critical  applications.  Furthermore,  it  has  resulted  in  high  material  qualification 
costs  and  overdesign  weight  penalties. 

Mechanistic  investigations  have  shown  that  several  processes  accompany  fracture 
of  composites.  From  a  physical  view  point,  some  of  them  are  related  to  microme¬ 
chanics,  while  others  to  macromechanics,  thus  involving  a  wide  range  of  scales.  A 
complete  characterization  in  space  and  time  of  these  processes  has  often  been  a 
formidable  experimental  task.  From  an  analytical  viewpoint  it  is  a  non-smooth 
problem,  involving  various  controlling  parameters  with  associated  spatial  and  tem¬ 
poral  scales. 

Among  the  various  mechanisms  with  a  significant  contribution  to  fracture  resis¬ 
tance  of  composites,  crack  bridging  by  intact  fibers  is  considered  to  be  very  impor¬ 
tant.  This  is  demonstrated  by  the  large  number  of  research  works  that  have  been 
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reported  in  the  literature.  Most  of  these  works  are  dedicated  to  deriving  analytical 
expressions  for  the  stresses  at  the  crack  tip  or  toughness  for  cracks  that  exhibit 
steady  or  non-steady  growth.  However,  experimental  studies  to  attest  these  models 
have  not  received  the  same  level  of  attention.  This  is  perhaps  due  to  the  highly 
complex  interactions  between  the  various  events  and  the  difficulty  to  isolate  their 
interactions. 

Another  mechanism  that  is  also  important  in  the  increased  fracture  resistance  of 
composites  is  an  interaction  of  the  crack  and  the  fibers  ahead  of  the  crack.  Studies 
on  the  interaction  effects  have  also  been  reported  albeit  to  a  lesser  extent.  The 
essentially  three  dimensional  problem  that  ensues  renders  an  analytical  solution 
very  difficult.  Thus  simplifying  assumptions  are  adopted  for  particular  cases  [1,2]. 
A  more  common  approach  has  been  to  study  the  interaction  between  a  crack  and 
inclusion  [3]  or  crack  near  or  on  a  bi-material  interface  [4  -  6]  using  two  dimensional 
configurations  to  derive  solutions  for  the  stress  and  strain  fields  as  well  as  for  the 
stress  intensity  factors.  These  solutions  offer  valuable  insights  in  the  interaction 
effects  but  their  applicability  to  a  fibrous  composite  system  may  be  limited  due  to 
various  length  scales  involved  in  the  vicinity  of  a  crack  tip  (Fig.  1). 

To  investigate  some  of  the  effects  of  fiber  spacing  and  fiber  type  on  strength 
and  fracture,  Botsis  and  co-workers  [7  -  10]  carried  out  research  on  composites  with 
controlled  fiber  spacing.  The  results  of  these  studies  have  shown  that  strength  and 
fracture  correlate  well  with  fiber  spacing.  This  approach  offers  certain  advantages. 
By  controlling  the  fiber  architecture,  the  properties  of  the  constituent  materials 
and  the  interface,  different  types  of  behaviours  can  be  simulated  experimentally. 
Thus  some  of  the  pertinent  parameters  (fiber  size,  spacing,  interface,  etc.,)  affect¬ 
ing  strength,  crack  growth  and  debond  evolution  can  be  identified  and  thoroughly 
characterized.  Results  from  such  studies  can  be  an  important  supplement  to  efforts 


2 


fx. 


Fig.  1.  Merging  micro-  and  macro-mechanics  g 


aimed  at  characterizing  the  response  of  composites  particularly  when  a  continuum 
approximation  of  the  reinforcement  is  not  realistic  (Fig.  1).  Such  cases  may  arise 
when  studying:  (a)  crack  initiation,  (b)  driving  forces  of  relatively  short  cracks, 
(c)  the  area  arovmd  the  tip  of  a  long  crack,  (d)  criteria  for  fracture,  etc.  More¬ 
over,  fiber  spacing  reflects  the  structure  of  the  composite  and  is  an  intermediate 
scale  since  it  is  larger  than  the  scales  involved  in  fiber  debonding,  fracture  and  pull 
out.  It  is  smaller,  however,  than  the  macro  scales  that  arise  in  macrocracking  and 
damage  zones.  Using  this  intermediate  length  scale  within  a  continuum  mechanics 
framework,  a  set  of  constitutive  equations  can  be  developed  that  incorporate  this 
parameter  which  reflects  the  reinforcement  better  than  volume  average  parameters. 

The  present  work  is  aimed  at  investigating  the  effects  of  bridging  fibers  and  fibers 
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in  the  path  of  the  crack  on  the  stress  field  and  stress  intensity  factor.  Emphasis  is 
placed  on  a  detailed  stress  analysis  within  a  zone  around  the  crack  tip  in  an  effort 
to  establish  its  boundaries,  (see  Fig.  1).  Outside  this  zone  homogenization  theories 
may  be  employed  that  consider  the  composite  as  an  orthotropic  material.  At  the 
boundary  of  this  zone,  appropriate  conditions  need  to  be  imposed  to  consistently 
“glue”  the  two  domains.  Emphasis  in  this  report  is  primarily  numerical  and  rather 
practical  in  nature.  We  incorporate  the  inherent  inhomogeneity  and  anisotropy  in 
our  framework  at  the  variational  level  and  look  at  material  mismatch,  interphase 
boundaries,  interaction  c#the  interface  boundaries  with  each  other  and  also  with 
high  stresses  at  the  crack  tip.  These  aforementioned  issues  render  the  ensuing 
problem  ‘non-smooth’  in  the  classical  sense  and  necessitate  the  solution  method 
to  be  numerical  in  nature  and  founded  on  sound  mathematical  footings.  We  have 
employed  the  finite  element  method  because  it  is  based  on  variational  principles  and 
is  a  powerful  numerical  technique  that  offers  a  unified  approach  to  such  multiphase 
problems. 

The  report  is  divided  into  three  main  sections  as  follows:  Section  2  presents 
the  variational  form  of  the  problem.  Section  3  develops  the  finite  element  formu¬ 
lation  of  the  proposed  model.  Section  4  gives  an  account  of  the  evaluation  of  the 
stress  intensity  factor.  A  detailed  set  of  numerical  simulations  employing  a  carefully 
designed  fiber  architecture  to  uniaxial  tension  is  presented  in  Section  5. 
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2.  The  Variational  Framework 


The  insitu  strength  of  a  fibrous  composite  is  not  only  a  function  of  the  prop¬ 
erties  of  each  of  its  constituent  materials,  it  also  depends  on  the  various  energy 
dissipation  mechanisms  that  lead  to  matrix  cracking  and/or  inelastic  deformations, 
fiber  debonding,  slipping,  and  their  ultimate  fracture.  Some  of  these  mechanisms 
such  as  elasto-plastic  deformations,  damage  and  cracking  in  the  matrix  and  necking 
and  fracture  in  the  bridging  fibers  can  successfully  be  addressed  at  the  constitutive 
level.  However,  fiber  debonding  and  pull-out  that  leads  to  strong  discontinuities  in 
the  displacement  field  is  related  to  the  interphase  kinematics,  and  an  appropriate 
way  to  handle  them  is  to  incorporate  them  in  the  variational  framework.  Our  ob¬ 
jective  here  is  to  develop  a  variationally  sound  approach  that  can  be  used  to  derive 
an  energy  based  criteria  for  crack  propagation  and  arrest  in  fibrous  composites  in 
the  presence  of  one  or  more  of  the  above  mentioned  mechanisms.  Furthermore,  an 
energy  based  approach  for  crack  propagation  and  arrest  can  readily  be  extended  to 
the  nonlinear  regime  of  material  and  geometric  response.  A  general  framework  for 
such  problems  is  provided  by  the  Hu-  Washizu  variational  principle  which  considers 
displacements,  strains  and  stresses  as  the  independent  field  variables. 

Let  C  be  a  bounded  open  set  with  piecewise  smooth  boundary  F,  Usd  > 
2  denotes  the  number  of  spatial  dimensions.  Let  F/  represent  the  fiber-matrix 
interface  with  and  be  the  displacement  fields  at  the  interface  of  the  fibers 
and  the  matrix.  We  can  write  an  energy  functional  for  the  hybrid  system  as 


n(e,<T,n, 


e  :  c  :  £dCl-\- 


X{u^-u^)dr-Ue^t{u) 


(1) 

where  u,  <t  and  e  are  the  unknown  displacement,  stress  and  strain  fields,  respec¬ 
tively,  and  A  is  a  Lagrange  multiplier.  Uexti'^)  represents  the  external  work  on  the 
composite  system.  The  first  term  in  (1)  is  the  stored  energy  function,  the  second 
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term  is  a  Lagrange  multiplier  (cr)  enforcement  of  the  strain-displacement  relations 
(e  =  V®w),  and  the  third  term  is  a  Lagrange  multiplier  enforcement  of  the  con¬ 
tinuity  of  displacement  field  at  the  fiber  matrix  interface  and  corresponds  to  the 
interface  energy.  The  Lagrange  multiplier  A  can  be  interpreted  as  a  force  quantity 
required  to  prevent  relative  displacement  at  the  interface  boundary.  In  the  present 
work,  we  assume  a  perfect  bond  between  the  fiber  and  the  matrix  and  thus  pre¬ 
clude  any  debonding  or  pull  out  of  fibers.  Consequently  the  third  term  becomes 
identically  zero  leading  to  the  following  3-field  energy  functional. 

n(e,  cr,M);=i  f  £:c:£dQ,+  f  cr  •  [ —  e]  dO  —  Ilext (■**)  (2) 

2  Jn  Jo, 

We  also  want  to  develop  a  general  framework  that  can  encompass  the  compress¬ 
ible  as  well  as  slightly  compressible  or  nearly-  incompressible  material  systems. 
For  example,  metal  matrix  composites  normally  show  a  compressible  elastic  and 
isochoric  (volume  preserving)  inelastic  response  while  polymer  based  composites 
show  a  volume  preserving  or  incompressible  behavior  even  at  small  strains.  It  is 
well  known  that  Galerkin  approximation  to  incompressible  or  nearly-incompressible 
elasticity  leads  to  eroneous  results  (Hughes  [ll]),e.g.,  mesh  locking  phenomenon, 
pressure  oscillations,  checker-board  modes,  to  name  a  few.  In  the  last  two  decades 
extensive  research  efforts  have  been  directed  to  this  issue  and  there  are  some  well 
established  techniques  that  can  quite  successfully  handle  the  incompressible  range 
of  the  response  spectrum  (  Hughes  [12,  13],  Nagtegaal  et.  al.  [14],  Simo  et.  al.  [15]) 

We  assume  the  strain  field  to  be  divided  into  deviatoric  and  volumetric  compo¬ 
nents 

e(u,  9?)  :  =  dev  [V^u]  H — 

'^'sd 

=  H — —  [^  —  divu]  1  (3) 

'^'sd 

where  u  :  Q  ^  is  the  displacement  field,  (p  :  Q  ^  TZ  is  the  spherical  part  of 
£  and  is  regarded  as  an  independent  field,  and  dev[-]  :=  (•)  —  ^tr(-)l  represents 
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the  deviatoric  part  of  the  indicated  argument.  The  stress  field  which  is  compatible 
with  the  assumed  strain  field  can  be  expressed  as 

(t{u,  (p,  p)  =  dev  [D£{u,  (/?)]  +  pi  (4) 

where  D  represents  the  material  constitutive  matrix  and  p:  ^  7^  is  an  indepen¬ 

dent  function  which  is  to  be  interpreted  as  the  hydrostatic  pressure.  By  substitut¬ 
ing  (3)  and  (4)  into  the  3-  field  functional  (2),  we  obtain  the  modified  functional 
U(u,ip,p) 

Il{u,(p,p)  =  ^  I  £{u,(p)  :  c  :  £{u,(p)dQ,  +  f  p{divu  —  (p)dQ  —  Ilext{u)  (5) 

^  Jq  Jq. 

We  denote  the  variations  in  volumetric  strain  p  and  pressure  p  hy  'll)  and  g, 
respectively,  all  of  which  are  only  required  to  be  square  integrable,  i.e.  L2(fi)  over 
the  domain.  As  no  boundary  conditions  are  applied  on  stresses  or  strains,  the 
admissible  spaces  of  solutions  and  corresponding  spaces  of  variations  coincide,  i.e., 
V  ;=  {p\p  €  [L2mh  and  P  :=  {p\p  G  [L2(f2)]}. 

The  boundary  of  domain  Q,  denoted  by  T,  is  assumed  to  be  piecewise  smooth. 
We  further  assume  that  T  is  decomposed  into  two  non-overlapping  subregions  T^ 
and  Th,  where  Tp  is  the  portion  of  the  boundary  with  prescribed  displacement 
field  g,  and  T/,.  is  the  portion  of  the  boundary  with  prescribed  tractions  h.  The 
admissible  spaces  for  the  displacement  field  and  its  variations  are 

<S  =  {n  =  (6) 

V  =  {a;  G  u):  (7) 

where  S  is  the  space  of  trial  displacements  and  V  is  the  associated  space  of  weighting 
functions,  respectively.  denotes  the  space  of  square-integrable  functions 

along  with  their  generalized  derivatives  defined  over  Q,  and  Hq{Q)  is  the  subset  of 
whose  members  satisfy  zero  essential  boundary  conditions. 


The  minimization  problem  consists  of  finding  the  stationary  point  of  (5)  and 
satisfying  the  prescribed  boundary  conditions  which  are  assumed  to  take  the  form 

u  =  g  Vcc  G  Fp  (8) 

(T-n  =  h  'ixEVh  (9) 

where  n  denotes  the  unit  outward  normal  to  boundary  F. 

The  stationarity  conditions  associated  with  (5)  yield  the  variational  form  of  the 
problem.  For  the  case  of  small  strains  and  static  loads,  and  assuming  the  fibrous 
composite  to  occupy  a  region  0,  in  the  formal  statement  of  the  variational  form 


may  be  expressed  as:  Given  /  :  Q,  — >•  g  :  F^  — h  :Th  — 

find  {u,  (p,p)  E  S  X  V  X  V  such  that  for  all  (a>,  il),q)EVx'DxV 

(V®w,  dev[D£(w,  (p)]  +  pi)  =  (oj,  /)  +  (w,  h)r  Vo;  G  V  (10) 

{q,divu  —  (p)  =  0  yqEL^{Q,)  (11) 

(■0,  -p  +  ^tr[£>e(«,  (p)])  =  0  VV’  G  L^(fl)  (12) 


where  (■,•)  denotes  the  I/2(fl)  inner  product,  and  (•,-)r  represents  the  T2(F)  inner 
product. 
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3.  Finite  Element  Framework. 


We  now  introduce  the  finite-dimensional  approximations  for  displacements, 
strains  and  stresses  as  C  »S,  C  V,T>^  C  V  and  C  V ,  respectively. 
The  discrete  displacement  field  and  its  variation  are  continuous  globally,  while  the 
discrete  stress  and  strain  fields  are  assumed  to  be  discontinuous  across  element 
boundaries.  As  a  result  of  this  discontinuous  approximations,  the  discrete  versions 
of  (11)  and  (12)  can  be  written  over  domain  of  an  individual  element  Q®  as 


(g®,divtt®  -  (p^)cie  =  0 

(13) 

(^®,  -p®  -1-  itr[r>er(w,<p)])ne  =  0 

(14) 

Consequently,  the  stress  and  strain  fields  can  be  eliminated  at  the  element  level, 
leading  to  a  modified  displacement  type  formulation.  With  this  end,  we  design 
discrete  projection  operator  for  the  strain  and  pressure  fields 

=  P'^(x)a^  (15) 

/  =  P'^{x)(3^  (16) 

where  P^(x)  =  {'yi(a;),  •  •  •  ,7„^(a;)},  is  a  set  of  prescribed  functions  and  o:®, 
yS®  e  7?.”'  are  local  element  parameters  corresponding  to  strains  and  stresses, 
respectively. 

We  define  the  discrete  divergence  operator  as 

div«®  =  (17) 

where  6®^  =  l^S®,  and  JB®  is  the  usual  strain-displacement  matrix  composed  of 
derivatives  of  shape  functions. 

We  also  define  a  mean  operator  over  element  domain  as 

if®  :=  [  P(x)P{xfdQ  (18) 
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and  assume  that  is  invertible  (this  condition  is  ensured  if  the  functions 
in  P(£c)  are  independent).  Since  any  e  V'^  is  discontinuous  across 
element  boundaries,  equation  (13)  implies 

(p^  =  P^{x)H^  ^  [  P{x)divu^df}  (19) 

for  e  =  1, 2,  •  •  • ,  Uei-  Inserting  (17)  into  (19)  we  finally  obtain 

=  P'^{x)H^~'  [  P{x)b^dQd^  (20) 

Thus,  the  effect  of  the  interpolations  (15)  and  (16)  is  to  define  a  modified  discrete 
divergence  operator  by  expression  (20),  and  is  expressed  as  =  divu^. 

Similarly  substituting  interpolations  (15)  and  (16)  in  (14),  we  obtain 

p^  =  P^{x)H^~'  f  P{x)hi[D£{u^,p^)]dn  (21) 

Also  by  inserting  (20)  into  the  definition  of  the  assumed  strain  field  (3),  we 
obtain  the  corresponding  discrete  version.  Now  substituting  the  discrete  version  of 
the  assumed  strain  field  together  with  expressions  (20)and  (21)  in  equation  (10) , 
we  get  a  displacement  version  of  the  discrete  variational  form.  It  thus  leads  to  an 
algebraic  system  of  equations  where  the  primary  variable  is  the  displacment  field 
alone. 
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4.  Evaluation  of  Stress  Intensity  Factor 

Let  be  the  strain  energy  in  equation  (2).  We  evaluate  this  strain  energy 

for  two  configurations  of  the  specimen,  i.e.,  for  two  different  crack  lengths  I2 
{I2  >  h).  From  the  two  crack  lengths  and  the  calculated  strain  energies,  the  energy 
release  rate  for  mode  1  crack  growth  is  calculated  as 

TTStrain  rrstrain 

where  and  are  the  calculated  strain  energies  at  crack  lengths  li  and 

I2,  respectively.  Using  linear  elastic  fracture  mechanics,  the  stress  intensity  factor 
for  mode  I  crack  is  then  computed  as 


/  E 

E  =  d  - - 2  ^  (plane  strain) 

(23) 

K  =  d EG  (plane  stress) 

(24) 

where  E  and  v  are  the  Young’s  modulus  and  Poisson’s  ratio  of  the  matrix  respec¬ 
tively. 
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5.  Numerical  Examples 


In  this  section,  we  present  various  numerical  experiments  on  2-D  plane  strain 
specimens  with  straight  fibers  aligned  along  the  loading  direction.  We  study  the  in¬ 
teraction  between  the  stress  field  around  the  crack  tip  and  the  inter-phase  stress  field 
that  exists  at  the  fiber-matrix  interface.  As  it  was  pointed  out  in  the  introduction, 
the  rows  of  fibers  are  replaced  by  layers  of  effective  material.  Furthermore,  these 
numerical  siumlations  also  represent  a  through-the-thickness  transverse  crack  in  a 
thick  multilayered  composite.  Figure  2  shows  the  specimen  geometry  with  a  trans¬ 
verse  notch  and  reinforcing  fibers.  Invoking  symmetry  in  the  geometry  and  loading, 
and  assuming  the  crack  to  grow  in  a  self-similar  fashion,  we  have  discretized  only 
half  of  the  specimen.  The  height  (H)  of  the  discrete  specimen  is  one-half  that  of  its 
width,  i.e.,  H  =  W/2.  The  crack  length  is  one  third  the  specimen  width,  I  =  W/Z. 
All  reinforcing  layers  have  a  constant  diameter  d  =  Z/30.  Figure  3  shows  a  typical 
finite  element  mesh  with  one  fiber  in  the  path  of  the  crack.  We  have  assumed  a 
perfect  bond  between  fiber  and  matrix  in  the  present  simulations.  The  ratio  of 
Young’s  Modulus  for  the  various  fibers  with  respect  to  the  matrix  Ef/Em  is  5,  15, 
25,  and  35.  The  Poisson’s  ratio  for  the  reinforcement  and  matrix  is  i^f  =  I'm  =  0-3. 

In  our  numerical  simulations,  we  are  controlling  the  fiber  spacing  in  the  trans¬ 
verse  direction  to  investigate  the  inter-relation  between  the  interaction  of  crack  and 
reinforcement  and  crack  bridging. 
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5.1  Shielding  Effect:  Proximity  of  fiber  and  reduction  in  SIF 

Reinforcing  fibers  in  the  path  of  a  crack  play  an  important  role  in  increasing 
the  strength  of  a  composite.  This  increase  in  strength  is  because  of  the  interaction 
between  the  high  stress  at  crack  tip  and  the  high  stress  gradients  at  fiber-matrix 
interface.  This  first  test  case  investigates  effect  of  proximity  of  a  fiber  ahead  of  the 
crack  on  the  reduction  in  SIF.  The  specimen  is  loaded  by  applying  2%  strain  at 
the  displacement  boundary.  In  this  test  case  we  have  considered  only  one  layer  of 
fibers  that  lie  in  the  path  of  the  crack.  The  distance  6  between  the  reinforcement 
and  the  crack  tip  is  varied  (in  multiples  of  the  diameter  of  the  fiber)  between  range 
d  <  6  <  lOd.  The  results  are  normalized  with  respect  to  the  SIF  obtained  from 
homogeneous  isotropic  specimen  subjected  to  same  loading  conditions.  As  can  be 
seen  in  Fig.  4,  bringing  the  fiber  close  to  the  crack  reduces  SIF  at  the  crack  tip. 
Furthermore,  a  stronger  fiber  results  in  a  greater  reduction  in  SIF  as  compared  to 
a  less  stiff  fiber  because  of  the  stiff  inter-phase  stress  boundary.  It  is  also  worth 
noticing  that  after  a  distance  of  lOd,  the  effect  of  the  reinforcement  on  Kf  has 
diminished  considerably  for  all  four  cases  of  Ef /Em- 
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Normalized  distance  in  terms  of  fiber  diameter  (d) 


Fig.  4.  Single  fiber  test.  Effect  of  fiber  proximity  on  reduction  in  K. 
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5.2  Effect  of  number  of  shielding  fibers  on  reduction  in  SIF 

This  test  case  investigates  the  effect  of  number  of  fibers  in  the  path  of  a  crack 
on  the  reduction  in  SIF.  The  results  of  the  simulations  are  shown  in  Fig.  5.  The 
first  reinforcing  fiber  is  one  diameter  away  from  crack  tip,  i.e.,  6  —  d.  The  center-to- 
center  distance  between  the  subsequent  fibers  is  2d,  i.e.,  A  =  2d.  Normalization  is 
performed  with  respect  to  the  homogeneous  isotropic  specimen  with  crack  of  same 
dimension.  The  data  in  Fig.  5  shows  that  the  first  fiber  has  the  greatest  effect  on 
Kf.  The  second  fiber  also  has  some  effect  but  the  third  and  fourth  fibers  have  a 
considerably  smaller  influence  on  Kf.  For  the  cases  of  Ef/Em  =  5  and  15  the  third 
and  fourth  fiber  have  practically  no  effect.  To  examine  any  fiber-fiber  interaction 
on  the  results  of  Fig.  5  we  performed  simulations  of  each  fiber  separately  at  their 
respective  locations.  The  result  of  the  simulations  showed  that,  for  the  parameters 
of  geometry  and  elastic  constants,  linear  superposition  can  be  used  for  the  overall 
effect  on  Kf.  To  investigate  the  effect  of  Ef/Em  on  Kf,  the  data  in  Fig.  5  are 
plotted  as  a  function  of  Ef/Em,  in  Fig.  6.  As  seen  in  Fig.  6,  the  trend  of  the 
dependence  is  very  similar  and  independent  of  the  number  of  fibers.  Apparently 
the  first  two  fibers  have  the  greatest  effect  on  the  total  stress  intensity  factor.  The 
addition  of  a  third  and  fourth  fibers  have  practically  no  effect. 
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5.3  Spatial  distribution  of  stress  field 

We  first  investigate  the  spatial  distribution  of  an  in  the  path  ahead  of  the  crack. 
Figure  7-a  shows  the  zoomed  view  of  crack  tip  singularity  in  the  homogeneous 
specimen.  Here  comparison  is  done  with  the  analytical  expressions.  We  know  that 
unlike  the  analytical  expression,  numerical  methods  always  give  a  finite  value  at 
crack  tip.  In  the  present  case,  the  numerical  method  does  capture  the  high  gradient 
at  crack  tip.  As  our  next  step,  we  investigate  the  reduction  in  the  intensity  of 
stress  at  crack  tip  due  to  reinforcing  fibers  in  the  path  of  ahead  of  crack.  Figure 
7-b  represents  this  reduction  in  crack  tip  stresses  due  to  two  shielding  fibers.  We 
can  clearly  see  the  drop  in  value  by  three  folds.  Furthermore,  we  can  see  the  stress 
gradients  across  the  width  of  the  fibers  as  well.  We  have  plotted  a  zoomed  view  of 
the  region  between  crack  tip  and  the  fibers  in  Fig.  7-c  which  shows  the  steep  stress 
gradients  at  the  fiber  matrix  interfaces.  We  have  also  plotted  ln(l/r)  singularity 
along  the  two  sides  of  adjacent  fibers.  Although  we  are  not  able  to  exactly  catch  the 
singularity,  we  certainly  capture  the  extremely  steep  gradients  in  the  stress  field. 
This  test  case  also  validates  the  robustness  of  the  formulation  being  used  here  for 
such  “mathematically  non-smooth”  problems. 
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Fig.  7a.  an  for  homogeneous  case. 


Fig.  7b.-  Crack-tip  and  fiber-matrix  singularities. 
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Fig.  7c.  Zoomed  view  of  fiber-matrix  singularities. 

To  gain  further  insight  on  the  stress  distribution  around  the  crack  tip  in  the 
presence  of  the  reinforcement  we  carried  out  two  sets  of  simulations.  In  the  first 
one  we  started  with  the  specimen  without  reinforcement.  Then  we  added  a  fiber  at 
distance  d  from  the  crack  tip  and  two  more  fibers  at  distances  d  from  the  first  one 
and  each  other.  We  have  not  considered  a  fourth  fiber  because  its  effect  on  Kf  was 
negligible  (Fig.  5).  We  have  loaded  the  specimen  by  applying  2%  strain. 
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Fig.  8a.  Effect  of  number  of  fibers  on  stress  field  -  Homogeneous  specimen 
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Fig.  8b.  Effect  of  number  of  fibers  on  stress  field  -  Fiber-matrix  Ef  jE.^ 


Figure  8  shows  selective  contours  of  the  an  stress  field,  which  are  markedly 
different  from  that  of  the  specimen  without  reinforcing  fibers.  We  can  also  see  the 
stress  contours  inside  the  fibers  and  extremely  high  stress  gradients  (sharp  boundary 
layers)  at  the  interface  of  the  two  materials.  The  iso-stress  lines  indicate  that  stress 
in  the  fibers  at  equal  distances  from  the  crack  path  depends  on  its  proximity  from 
the  crack  tip;  the  closest  fiber  being  under  the  most  stress.  In  this  simulation 
Ef/Em  ratio  is  kept  constant  and  equal  to  15. 
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In  the  second  simulation,  we  examined  the  effects  of  Ef/Em,  on  the  stress  fields 
CTii  around  the  crack  tip  in  the  presence  of  three  fibers.  Figure  9  presents  the  spatial 
stress  distribution  an  for  Ef/Em  equal  to  5  and  25,  while  keeping  the  number  of 
fibers  fixed  equal  to  3.  As  can  be  seen,  the  inter-phase  stress  boundary  layers  of 
stronger  fibers  e.g..  Fig.  9-b  have  a  larger  zone  of  influence  as  compared  to  that  of 
the  weaker  fibers,  i.e..  Fig.  9-a.  Furthermore,  interaction  of  the  inter-phase  stress 
boundaries  can  be  seen  in  these  figures. 
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5.4  Effect  of  a  bridging  fiber 

In  this  numerical  simulation  we  introduce  an  intact  fiber  in  the  bridging  zone 
in  addition  to  one  reinforcing  fiber  in  the  path  ahead  of  the  crack  tip  (see  Fig.  10). 
The  motivation  for  this  case  study  originates  from  the  need  to  examine  the  relative 
effects  of  a  bridging  fiber  and  a  fiber  ahead  of  the  crack  tip  on  the  stress  intensity 
factor.  To  examine  the  contribution  of  a  bridging  fiber  we  vary  the  distance  A 
between  the  fiber  and  the  crack  tip  in  terms  of  multiples  of  fiber  diameter  d  with 
d  <  A  <  6d.  A  schematic  of  the  fiber  location  is  shown  in  Figure  10.  The  results 
of  the  simulations  are  shown  in  Fig.  11.  Note  that  the  reduction  in  normalized  SIF 
is  inversely  proportional  to  fiber  proximity  in  the  wake  of  the  crack.  Normalization 
here  is  done  with  respect  to  homogeneous  isotropic  specimen  with  a  crack  of  same 
dimension.  Comparing  Fig.  4  and  Fig.  11  for  Ef/Em=  5  and  S  =  d,  we  see 
that  a  fiber  in  the  wake  of  crack  reduces  the  normalized  SIF  to  3%  of  its  values  as 
compared  to  the  one  ahead  of  crack  tip  which  reduces  it  down  to  approximately 
87%.  Consequently,  for  given  Ef /Em.  ratio  and  given  8  from  the  crack  tip,  a  bridging 
fiber  is  much  more  effective  as  compared  to  its  counterpart  ahead  of  the  crack  tip 
at  the  same  distance. 
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Fig.  10.  Schematic  diagram  of  the  problem. 

The  normalized  force  in  the  bridging  fiber  (at  the  crack  face)  as  a  function  of 
EffEm  ratio  for  various  A  values  is  plotted  in  Fig.  12.  Here  normalization  in  each 
case  of  A  is  done  by  dividing  the  force  for  various  Ef  / Em.  ratios  by  their  correspond¬ 
ing  value  of  force  for  Ef/Em.  =  5.  It  can  be  seen  that  for  a  given  fiber  crack-tip 
distance  A,  the  stronger  fibers  carry  more  stress  and  shift  the  concentration  of  stress 
away  from  crack  tip,  thus  resulting  in  considerable  reduction  in  SIF. 
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Normalized 


5.5  Effects  of  the  number  of  fibers  in  the  bridging  zone 

To  examine  the  reduction  in  SIF  due  to  the  number  of  fibers  in  the  bridging  zone 
and  stress  level  along  the  crack  path,  we  carried  out  simulations  with  three  fibers 
ahead  of  the  crack  tip  and  adding  one  fiber  at  a  time  in  the  bridging  zone.  We  have 
taken  only  three  fibers  along  the  crack  path  because  the  effects  of  additional  fibers 
reduces  the  SIF  insignificantly  (Fig.  5).  Loading  is  the  same  as  in  the  previous 
tests.  The  distance  of  the  first  fiber,  ahead  of  the  crack  tip  is  one  diameter,  i.e., 
6  =  d,  and  the  surface-to-surface  interfiber  distance  is  two  diameters,  i.e.,  A  =  2d. 
The  distance  of  the  first  bridging  fiber  from  crack  tip  is  Id,  i.e.,  A  =  d.  The  fibers 
in  the  bridging  zone  are  added  at  distances  of  2d  (center-to-center). 


Number  of  brigding  fibers 

Fig.  13.  Effect  of  number  of  bridging  fibers  on  reduction  in  K. 


Figure  13  shows  the  evolution  of  the  normalized  SIF  as  a  function  of  the  number 
of  fibers  in  the  bridging  zone  for  the  case  when  Ef/Em  =  25.  The  numerical 
data  clearly  demonstrates  that  the  additional  fibers  further  decrease  the  SIF  at  the 
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crack  tip.  However,  the  contribution  of  the  first  fibers  is  noticeably  larger  than 
that  of  other  fibers.  Moreover,  the  infiuence  starts  diminishing  substantially  after 
the  fifth  fiber.  The  effects  of  the  ratio  Ef/Em  on  the  SIF  are  displayed  in  Fig. 
14.  It  is  interesting  to  notice  that  the  effects  of  Ef/Em  diminish  rapidly  and  after 
Ef/Em.  -  15  any  further  increase  in  Ef/Em  does  not  infiuence  the  SIF. 


Fig.  14.  Effect  of  relative  strength  of  the  fiber-matrix  system  on  reduction  in  K. 
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In  this  work  we  have  also  studied  the  normal  stress  (along  line  x  =  0,  see  e.g. 
Fig.  10)  ahead  of  the  crack  and  in  the  bridging  zone.  The  data  of  these  simulations 
are  shown  in  Fig.  15  for  Ef/Em  —  15  (a)  and  Ef/Em  —  25  (b).  Notice  here  that 
the  stress  in  the  fibers  ahead  of  the  crack  tip  is  approximately  the  same  in  all  three 
fibers,  and  about  the  same  level  as  that  given  by  the  rule  of  mixtures.  As  for  the 
stresses  in  the  fibers  in  the  bridging  zone,  there  seems  to  be  an  small  increase  of  the 
stress  except  for  the  fiber  located  the  furthest  from  the  crack  tip.  This  fiber  suffers 
the  largest  stress  but  contributes  the  least  on  the  SIF. 

Figure  16  shows  the  evolution  of  an  on  the  fibers  in  the  bridging  zone.  It  is 
interesting  to  notice  that  a  cut-off  length  where  the  stress  on  the  fibers  equals  that 
given  by  the  rule  of  mixtures  decays  very  quickly  as  we  approach  the  crack  tip. 
Furthermore,  it  can  be  seen  that  the  computed  stress  field  does  not  exhibit  any 
spurious  oscillations  or  checker  board  modes  even  under  extremely  high  gradients. 
This  simulation  makes  it  precise  that  in  order  to  study  fibrous  composites,  the 
numerical  method  should  be  based  on  sound  variational  foundations  because  we 
inevitably  have  to  use  energy  based  criteria  for  crack  propagation  in  these  materials. 
Furthermore,  it  should  also  be  able  to  capture  the  stresses  as  accurately  as  possible 
and  should  have  the  capacity  to  model  steep  gradients  in  stress  fields. 
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Fig.  16.  3-D  projection  of  an  field. 


32 


5.6  Effect  of  uniform  thermal  stresses 

Many  practical  applications  of  fibrous  composites  are  under  extreme  thermal  en¬ 
vironments.  Consequently,  it  seems  appropriate  to  investigate  the  effects  of  thermal 
stresses  on  the  stress  field  at  crack  tip  in  the  matrix. 

Temperature  effects  can  be  incorporated  in  the  variational  framework  if  the 
distribution  of  change  in  temperature  AT(x)  is  known. 

AT{x)=T{x)-To{x) 

where  T{x)  is  the  current  spatial  distribution  of  temperature  and  To{x)  is  the 
reference  temperature  distribution  for  zero  thermal  stresses  in  the  composite.  Let 
be  the  pointwise  thermal  gradient  vector  defined  as 

=  K{axAT  ayAT  0)^ 

where  k  is  1  for  plane  stress  and  (1  -I-  u)  for  plane  strain,  ax  and  ay  are  the 
coefficients  of  thermal  expansion  in  x  and  y  directions,  respectively.  Stresses  are 
related  to  thermal  strains  via  constitutive  relation  <t  =  c  :  {e  —  e***).  The  effect 
of  temperature  can  be  accounted  for  at  the  variational  level  by  replacing  the  first 
integral  in  equation  (2)  with  a  modified  strain  energy  term  written  as 

nstrain  =  ^  /  ■  C  :  {s  - 

^  Jn 

=  i  [  [{e:c:e)  +  (e***  :  c  :  -  2  (e  :  c  :  e**')]  dCl  (25) 

2  Jq 

where  the  first  term  in  (25)  is  the  usual  strain  energy  and  identical  to  the  first 
integral  in  (2)  while  the  remaining  two  terms  incorporate  the  temperature  effects. 

The  present  simulation  considers  the  thermal  stresses  that  arise  because  of  a 
change  in  the  temperature  of  the  fiber  matrix  system.  In  these  simulations  we 
load  the  specimen  by  applying  a  specified  force  at  the  boundary.  The  ratio  of 
thermal  coefficient  for  fiber  and  matrix  a/ /urn  =  10“^  and  is  assumed  to  be  constant 
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over  the  range  of  temperature  considered.  The  specimen  consists  of  one  bridging 
fiber  a.t  6  =  2d  and  no  fibers  in  the  path  ahead  of  the  crack.  Figure  17  shows 
a  linear  reduction  in  SIF  for  various  Ef/Em  ratios.  We  have  also  plotted  the 
normalized  force  at  the  fiber  tip  (Fig.  18)  where  normalization  is  done  with  respect 
to  zero  temperature  case.  It  can  be  observed  that  there  is  a  corresponding  linear 
increase  in  force  in  the  fiber.  This  simulation  points  out  that  a  variable  coefficient  of 
thermal  expansion  can  play  an  important  role  in  the  analysis  of  bimaterial  specimen. 
Consequently,  a  fully  coupled  thermo-mechanical  theory  that  would  evaluate  the 
spatial  distribution  of  thermal  gradients  and  consequently  effect  the  resultant  stress 
field  pointwise  would  be  more  appropriate  for  the  analysis  of  fibrous  composite 
systems. 
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Normalized  force  in  fibers 


Fig.  18.  Normalized  reaction  in  the  bridging  fiber  vs.  temperature. 


References 


[1]  H.  Gao,  J.R.  Rice,  “A  first-order  perturbation  analysis  of  crack  trapping  by 
arrays  of  obstacles”,  ASME  Journal  of  Applied  Mechanics,  56,  828,  1989 

[2]  A.F.  Bower,  M.  Ortiz,  “A  three-dimensional  analysis  of  crack  trapping  and 
bridging  by  tough  particles” ,  Journal  of  the  Mechanics  and  Physics  of  Solids, 
39,  815,  1991 

[3]  A. A.  Rubinstein,  “Macrocrack-microdefect  interaction’,  ASME  Journal  of  Ap¬ 
plied  Mechanics,  53,  505,  1986 

[4]  C.  Atkinson,  “On  the  stress  Intensity  Factors  Associated  with  Cracks  Inter¬ 
acting  with  an  Interface  Between  wo  Elastic  Media”,  International  Journal  opf 
Engineering  Science,  13,  487-504,  1975 

[5]  F.  Erdogan,  G.D.  Gupta  and  T.S.  Cook,  “Numerical  solutions  of  Singular  Inter- 
gral  Equation”,  1,  ed.  G.C.  Sih,  Noordhoff  International  Publishing,  1053-1067, 
1973 

[6]  A.  Romeo,  R.  Ballarini,  “A  Crcack  Very  Close  to  a  Bimaterial  Interface”,  AS- 
MEJournal  of  Applied  Mechanics,  62,  614  -619,  1995 

[7]  J.  Botsis,  C.  Beldica  and  D.  Zhao,  International  Journal  of  Fracture,  63,  375, 

1995 

[8]  D.  Zhao,  J.  Botsis,  “Experimental  and  Numerical  Studies  in  Model  composites 
Part  I  -  Experimental  Results”,  International  Journal  of  Fracture,  82,  153-174, 

1996 

[9]  C.  Beldica,  J.  Botsis,  “Experimental  and  Numerical  Studies  in  Model  compos¬ 
ites  Part  II  -  Numerical  Results”,  Journal  of  Fracture,  82,  175-192,  1996 

[10]  J.  Botsis,  D.  Zhao,  ’Fatigue  Fracture  Processes  in  a  Model  Composite’,  Com¬ 
posites,  Part  A,  28A,  657-666,  1997 

[11]  T.J.R.  Hughes,  The  Finite  Element  Method:  Linear  Static  and  Dynamic  Finite 
Element  Analysis,  Prentice-Hall,  Englewoods  Cliffs,  NJ,  1987. 

[12]  T.J.R.  Hughes,  “Generalization  of  selective  integration  procedures  to 
anisotropic  and  nonlinear  media” ,  International  Journal  for  Numerical  Methods 
in  Engineering,  15,  1413-1418,  1980 

[13]  T.J.R.  Hughes,  “Recent  developments  in  computer  methods  for  structural  anal¬ 
ysis”,  Nuclear  Engineering  and  Design,  57  No.  2,  427-439,  1980 


36 


[14]  J.C.  Nagtegaal,  D.M.  Parks,  and  J.R.  Rice,  “On  numerically  accurate  finite 
element”.  Computer  Methods  in  Applied  Mechanics  and  Engineering  ,4,  153- 
178,  1974 

[15]  J.C.  Simo,  R.L.  Taylor,  and  K.S.  Pister,  “Variational  and  projection  meth¬ 
ods  for  the  volume  constraint  in  finte  deformation  elasto-plasticty” ,  Computer 
Methods  in  Applied  Mechanics  and  Engineering,  51,  177-208,  1985 


37 


